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ABSTRACT 

The evolution of a system consisting of a protoplanetary disc with two embedded 
Jupiter sized planets is studied numerically. The disc is assumed to be flat and non- 
self gravitating, which is modeled by the planar (two-dimensional) Navier-Stokes equa- 
tions. The mutual gravitational interaction of the planets and the star, and the gravi- 
tational torques of the disc acting on the planets and the central star are included. The 
planets have an initial mass of one Jupiter mass Mj^p each and the radial distances 
from the star are one and two senii-niajor axis of Jupiter, respectively. 

During the evolution both planets increase their mass due to accretion of gas from 
the disc; after about 2500 orbital periods of the inner planet it has reached a mass of 
2.3 and the outer planet of 3.2 Mjup- The net gravitational torques exerted by the 
disc on the planets result in an inward migration of the outer planet on time-scales 
comparable to the viscous evolution time of the disc, while the semi- major axis of 
the inner planet remains constant. When the distance of close approach eventually 
becomes smaller than their mutual Hill radius the eccentricities increase strongly and 
the system may turn unstable. 

The implications for the origin of the solar system and the extrasolar planets are 
discussed. 
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It is generally assumed that planetary systems form in a dif- 
ferentially rotating gaseous disc. In the late stages of their 
formation the protoplanets are still embedded in the proto- 
stellar disc and their orbital evolution is coupled to that of 
the disc. Gravitational interaction between the planets and 
the gaseous disc has basically two effects: 

a) The torques by the planets acting on the disc tend to push 
away material from the orbital radius of the planet, and for 
sufficiently massive planets a gap is formed in the disc (Pa- 
paloizou & Lin 1984; Lin & Papaloizou 1993). The dynami- 
cal process of gap formation has been studied through time 
dependent hydrodynamical simulations for planets on circu- 
lar orbits by Bryden et al. (1999) and Kley (1999, henceforth 
paper I). The results indicate that even after the formation 
of a gap, the planet may still accrete material from the disc 
and reach about 10 Jupiter masses {Mjup)- For very low 
disc viscosity and larger planetary masses the mass accu- 
mulation finally terminates (Bryden et al. 1999). 

b) Additionally, the gravitational force exerted by the disc 
alters the orbital parameter (semi-major axis, eccentricity) 
of the the planet. Here, these forces typically induce some in- 
ward migration of the planet (Goldreich & Tremaine 1980) 



which is coupled to the viscous evolution of the disc (Lin 
& Papaloizou 1986) . Hence, the present location of the ob- 
served planets (solar and extrasolar) may not be identical 
with the position at which they formed. 

In particular, the migration scenario applies to some of 
the extra-solar planets (for a summary of their properties see 
Marcy, Cochran & Mayor 1999), the 51 Peg-type planets. 
They all have masses of the order Mjup, and orbit their 
central stars very closely, having orbital periods of only a few 
days. As massive planets, according to standard theory, have 
formed at a few AU distance from their stars, these planets 
must have migrated to their present position. The inward 
migration was eventually halted by tidal interaction with the 
star or through interaction with the stellar magnetosphere 
(Lin, Bodenheimer & Richardson 1996). The only extrasolar 
planetary system known so far {v And) consists of one planet 
at 0.059 AU on a nearly circular orbit and two planets at .83 
and 2.5 AU having larger eccentricities (.18 and .41) (Butler 
et al. 1999). 

In case of the solar system the question, what prevented 
any further inward migration of Jupiter, arises. As the net 
tidal torque on the planet is a delicate balance between the 
torque of the material located outside of the planet and the 
material inside (eg. Ward 1997), any perturbation in the 
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density distribution may change this balance. In this letter 
we consider the effect that an additional planet in the disc 
has on the migration rate. 

We present the results of numerical calculations of a 
thin, non-self gravitating, viscous disc with two embedded 
protoplanets. Initially the planets with one Mjup each are 
on circular orbits at a = lajup and 2 ajup, respectively. In 
contrast to the existing time-dependent models (Bryden et 
al. 1999; paper I) we take into account the back-reaction of 
the gravitational force of the disc on the orbital elements 
of the planets and star. The models are run for about 3000 
orbital periods of the inner planet corresponding to 32,000 
years. In Section 2 a description of the model is given, the 
results are presented in Section 3 and our conclusions are 
given in Section 4. 



2 THE MODEL 

We consider a non-self-gravitating, thin accretion disc model 
for the protoplanetary disc located in the z = plane and 
rotating around the 2-axis. Its evolution is described by the 
two dimensional (r — ip) Navier-Stokes equations, which are 
given in detail in Kley (1999, paper I). The motion of the disc 
takes place in the gravitational field of the central star with 
mass Ms and the two embedded protoplanets with masses 
m\ and 7712. In contrast to paper I we use here a non-rotating 
frame as both planets have to be moved through the grid. 
The gravitational potential is then given by 
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where G is the gravitational constant and r^, ri, and r2 
are the radius vectors to the star and two planets, respec- 
tively. The quantities si and S2 are smoothing lengths which 
are 1/5 of the corresponding sizes of the Roche-lobes. This 
smoothening of the potential allows the motion of the plan- 
ets through the computational grid. 

The motion of the star and the planets is determined 
firstly by their mutual gravitational interaction and secondly 
by the gravitational forces exerted on them by the disc. The 
acceleration of the star as is given for example by 
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where the integration is over the whole disc surface, and E 
denotes the surface density of the disc. The expressions for 
the two planets follow similarly. We work here in an accel- 
erated coordinate frame where the origin is located in the 
centre of the (moving) star. Thus, in addition to the gravita- 
tional potential (|l|) the disc and planets feel the additional 
acceleration —as. 

The mass accreted from the disc by the planets (see 
below) has some net angular momentum which in principle 
changes also the orbital parameter of the planets. However, 
this contribution is typically about an order of magnitude 
smaller than the tidal torque (Lin et al. 1999) and is ne- 
glected here. 

As the details of the origin and magnitude of the vis- 
cosity in discs is still uncertain we assume a Reynolds-stress 



formulation (paper I) with a constant kinematic viscosity. 
The temperature distribution of the disc is fixed through- 
out the computation and is given by the assumption of a 
constant ratio of the vertical thickness H and the radius r. 
Hence, the fixed temperature profile is given by T{r) oc r~^. 
We assume H/r = 0.05, which is equivalent to a fixed Mach 
number of 20. 

For numerical convenience we introduce dimensionless 
units, in which the unit of length, Ro, is given by the initial 
distance of the first planet to the star Rq = ri(t = 0) = 
lajup. The unit of time is obtained from the (initial) orbital 
angular frequency fii of the first planet i.e. the orbital period 
of the planet 1 is given by 



Pi = 2-nto. 



(3) 



The evolutionary time of the results of the calculations as 
given below will usually be stated in units of Pi. The unit 
of velocity is then given by vq = Ro/to- The unit of the 
kinematic viscosity coefficient is given by vo = Ro^o- Here 
we take a typical dimensionless value of 10~^ corresponding 
to an effective a of 4 x 10~^. 



2.1 The numerical method in brief 

The normalized equations of motion are solved using an Eu- 
lerian finite difference scheme, where the computational do- 
main [r,nin,rmax] X [ifimin , ifimax] is Subdivided into A'"r X N^ 
grid cells. For the typical runs we use Nr = 128, N,p = 128, 
where the azimuthal spacing is equidistant, and the radial 
points have a closer spacing near the inner radius. The nu- 
merical method is based on a spatially second order accu- 
rate upwind scheme (monotonic transport), which uses a 
formally first order time-stepping procedure. The method- 
ology of the finite difference method for disc calculations is 
outlined in Kley (1989) and paper I. 

The N-body module of the programme uses a forth or- 
der Runge-Kutta method for the integration of the equations 
of motion. It has been tested for long term integrations us- 
ing the onset of instability in the 3-body problem consisting 
of two closely spaced planets orbiting a star as described by 
Gladman (1993). For the initial parameter used here, the 
error in the total energy after 1.2 x 10"" orbits (integration 
over lO^yrs) is less 2 x 10~^. 



2.2 Boundary and initial conditions 

To cover the range of influence of the planet on the disc 
fully, we typically choose for the radial extent (in dimen- 
sionless units, where planet 1 is located initially at r = 1) 
rmin ~ 0.25, rmax = 4.0. The azimuthal range covers a com- 
plete ring (fmin = 0.0, (fmax = ^Tv usiug pcriodic boundary 
conditions. To test the accuracy of the migration, a com- 
parison calculation with r,nax = 8.0 and higher resolution 
Nr = 256, Nip — 256 was also performed. 

The outer radial boundary is closed to any mass flow 
«(fmax) = 0, while at the inner boundary mass outflow is 
allowed, emulating accretion onto the central star. At the 
inner and outer boundary the angular velocity is set to the 
value of the unperturbed Keplerian disc. Initially, the matter 
in the domain is distributed axially symmetric with a radial 
surface density proflle E ex r~^". 
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Figure 1. The azimuthally averaged surface density for different 
times (stated in units of Pi). The solid Une refers to the initial 
density distribution. The density inside planet 1 is reduced be- 
cause of mass leaving through r^in . 



Figure 2. Grey-scale plot of the surface density of the at four 
different times in an 93 — r diagram. The small ellipses indicate 
the positions of the planets and sizes of their Roche lobes. 



Two planets, each with an initial mass of IMjup, are 
located at ri = 1.0, (pi = tt and r2 = 2.0, 1^2 = 0. Thus, they 
are not only spaced in radius but are positioned (in ip) in 
opposition to each other to minimize the initial disturbance. 
The radial velocity v is set to zero, and the angular velocity 
is set to the Keplerian value of the unperturbed disc. 

Around the planets we then introduce an initial den- 
sity reduction whose approximate extension is obtained from 
their masses and the magnitude of the viscosity. This initial 
lowering of the density is assumed to be axisymmetric; the 
radial profile E(r) of the initial distribution is displayed in 
Fig. 1 (solid line). The total mass in the disc depends on the 
physical extent of the computational domain. Here we as- 
sume a total disc mass within rmin = 0.25 and Vmax ~ 4.00 
of 0.01 Mq. The starting model is then evolved in time and 
the accretion rates onto the planets is monitored, where a 
given fraction of the mass inside the Roche-lobe of the planet 
is assumed to accrete onto the planet at each time step and 
is taken out of the computation and added to masses of 
the planets. The Courant condition yields a time step of 
6.810-"* Pi. 



3 RESULTS 

Starting from the initial configuration (Fig. 1) the planets 
exert torques on the adjacent disc material which tend to 
push mass away from the location of the planets. At the 
same time the planets continuously accrete mass from its 
surroundings, and the mass contained initially between the 
two planets is quickly accreted by the planets and added 
to their mass. Finally one large gap remains in the region 
between r = 1 and r = 2 (Fig. 1). 

Similarly to the one planet calculations as described in 
Bryden et al. (1999) and in paper I, each of the two plan- 



ets creates a spiral wave pattern (trailing shocks) in the 
density of the disc. In case of one disturber on a circular 
orbit the pattern is stationary in the frame co-rotating with 
the planet. The presence of a second planets makes the spi- 
rals non-stationary as is seen in the snapshots after 50, 100, 
250 and 500 orbits of the inner planet that are displayed 
in figure 2. Near the outer boundary at r = 4 the reflec- 
tion of the spiral waves are visible. Using the higher resolu- 
tion model with rmax = 8.0 (section /refbounds) we tested 
whether the numerical resolution or the wave reflection at 
rmax has any influence on the calculation of the net torques 
acting on the planet and the accretion rates onto the planet. 
Due to limiting computational resources the higher resolu- 
tion model was run only for a few hundred orbital periods 
and the largest difference (2.5%) occurred in the mass ms 
of the outer planet. The difference in radial position (migra- 
tion) is less than 1%. We may conclude that our resolution 
was chosen sufficiently fine and that the reflections at the 
outer boundary rmax = 4 do not change our conclusions 
significantly. 

In previous calculations (paper I) the equilibrium mass 
accretion rate from the outer part of the disc onto a one 
Jupiter mass planet for the same viscosity {u = 10"^) and 
distance from the star was found to be 4.35 x lO^'^Mjup/yr 
for a fully developed gap. Here the accretion rate onto 
the planets is much higher in the beginning (~ 5 x 
10~*Mjup/yr) as the initial gap was not as cleared. Thus, 
during this gap clearing process, the masses of the individual 
planets grows rapidly at the onset of calculations (Fig. 3). 
At i ~ 250 the mass within the the gap has been exhausted 
(see also Fig. 1) and the accretion rates M on the plan- 
ets lower dramatically. At later times after several hundred 
orbits (Pi) they settle to nearly constant values of about 
10~^Mjup/yr for the outer planet, and 2.2 x 10~'^ Mjup/yr 
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Figure 3. The increase of the masses of the planets in units of 



M 



Jup- 



1.6 



Planet>x^|^ 








^^^ 


'**««»*»^ 


'™''*"^*ei««iS§a&a!i 


*^ 
















Planet 1 


gllllllll 


llllllllj 


1 



1000 1500 

Time 



Figure 4. The evolution of the semi-major axis of the planets. 



for the inner planet (from Fig. 3). Since the mass inside of 
planet 1 has left the computational domain and the initial 
mass between the two planets has been consumed by the 
two planets, this mass accretion rate onto planet 1 for later 
times represents the mass flow of material coming from radii 
larger than r2 (beyond the outer planet). It is the material 
which has been flowing accross the gap of the outer planet. 
Previously, this mass flow across a gap has been calculated 
to be about one seventh of the mass accretion rate onto a 
planet (paper I) and the present results are entirely consis- 
tent with that estimate. 

The gravitational torques exerted by the disc lead to 
an additional acceleration of the planets resulting in an ex- 
pression similar to the acceleration of the star (Eq. Q). For 
one individual planet this force typically results in an in- 
ward migration of the planet on timescales of the order of 
the viscous time of the disc (Lin & Papaloizou 1986) . Here 
this inward migration is seen clearly for the outer planet in 
Fig. 4 where the time evolution of the semi-major axis of 
the two planets is plotted. 

The inner planet on the other hand initially, for t < 200, 
moves slightly inwards but then the semi-major axis in- 
creases and, showing no clear sign of migration anymore, 
settles to a constant mean value of 1.02. However, the de- 
crease of the semi-major axis of the outer planet reduces the 
orbital distance between the two planets. From three body 
simulations (a star with two planets) and analytical consid- 
erations (Gladman 1993; Chambers, Wetherill & Boss 1996) 
it is known that when the orbital distance of two planets lies 
below the critical value of Acr ~ 2\/^Rh, where 



Rh = 



mi + 7712 \ ^^^ ai + a2 



/mi -f 7712 \ 

V 3M, J 



(4) 



is the mutual Hill radius of the planet, the orbits of the 
planets are not stable anymore. In the calculations this ef- 
fect is seen in the strong increase of the eccentricity of the 
inner planet. At i = 2500 its eccentricity has grown to about 
ei — 0.1, while the eccentricity of the outer planet remains 
approximately constant at a level of 62 — 0.03. 

We should remark here that in the pure 3-body prob- 
lem without any disc and the same initial conditions (ri — 
1.0,7711 — 1.0; r2 — 2.0,7712 = 1-0) for the three objects, the 



semi- major axis of the planets stay constant as this system is 
deflnitely Hill stable (Gladman 1993). However, if one takes 
as initial conditions for the pure 3-body system the param- 
eters for the planets as obtained from the disc evolution at 
t = 2500 (ri = 1.0, mi = 2.3; r2 = 1.5, m2 = 3.2) then the 
evolution becomes chaotic on timescales of hundreds of or- 
bits and the eccentricity grows up to e = 0.6 for both planets 
within 4000 orbits. 



4 CONCLUSIONS 

We have presented calculations of the long term evolution of 
two embedded planets in a protoplanetary disc that covers 
several thousand orbital periods. The planets were initially 
located at one and two ajup from the central star with initial 
masses of IMjup each. The gravitational interaction with 
the gaseous disc having a total mass of O.OIM© leads to an 
inward migration of the outer planet while the semi-major 
axis of the inner planet remains approximately constant and 
even slightly increased. At the same time, the ongoing ac- 
cretion onto the planets increases their masses continuously 
until at the end of the simulation (at t — 2500) the outer 
planet has reached a mass of about 3.2Mj„p and the inner 
planet of about 2.3Mjup- 

This increase in mass and the decreasing distance be- 
tween them renders the orbits eventually unstable resulting 
in a strong increase of the eccentricities on timescales of a 
few hundred orbits. 

From the computations we may draw three major con- 
clusions: 

1) The inward migration of planets immersed in an accretion 
disc may be halted by the presence of additional protoplan- 
ets located for example at larger radii. They disturb the 
density distribution significantly which in turn reduces the 
net gravitational torque acting on the inner planet. Thus, 
the migration of the inner planet is halted, and its semi- 
major axis remains nearly constant. 

2) When disc depletion occurs sufficiently rapid to prevent 
a large inward migration of the outer planet(s), a planetary 
system with massive planets at a distance of several au re- 
mains. This scenario may explain why for example in the 
solar system the outer planets (in particular Jupiter) have 
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not migrated any closer towards the sun. 
3) If the initial mass contained in the disc is sufficiently 
large then the inward migration of the outer planet (s) will 
continue until some of them reach very close spatial separa- 
tions. This will lead to unstable orbits resulting in a strong 
increase of the eccentricities. Orbits may cross and plan- 
ets may then be driven either to highly eccentric orbits or 
may leave the system all together (see eg. Weidenschilling & 
Marzari 1996). This may then explain the high eccentricities 
in some of the observed extrasolar planets in particular the 
planetary system of v Andromedae. 

As planetary systems containing more than two planets 
display different stability properties (Chambers et al. 1996) 
it will be interesting to study the evolution of multiple em- 
bryos in the protoplanetary nebula. 
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